%%
N = 5;
H = 2^N-1;
Z = zeros(N, H);
for h = 1:H
    bcode = dec2bin(h)-'0';
    Z(N-length(bcode)+1:N, h) = bcode';
end

K = rand(H,1) > 0.5;
KK = diag(K);

D = 1;
mu = rand(1, D);
SIGMA = rand(D,D);
SIGMA = (SIGMA*SIGMA') + eye(D);

%%
nsamples = 10000;
XX = zeros(N, nsamples);
for i = 1:nsamples
    V = mvnrnd(mu, SIGMA, H)';
    X = Z*KK*V';
    XX(:,i) = X;
end
covX = cov(XX');

%%
covtimes = covX./SIGMA;
rndcovtimes = round(covtimes);
assert(norm(rndcovtimes-rndcovtimes', 'fro') == 0);
datastats(covtimes(:) - rndcovtimes(:))
covpat = zeros((N^2+N)/2, H);
n = 1;
for i = 1:N
    for j = i:N
        covpat(n,:) = Z(i,:) & Z(j,:);
        n = n+1;
    end
end
covtimesv = nonzeros(tril(rndcovtimes));
assert(~any(covpat*K - covtimesv));
Kp = covpat\covtimesv;
rndKp = round(Kp);

%%
figure, bar([K, rndKp]);
figure, bar(covpat*(rndKp-K));